Quantum master equation scheme of time-dependent density functional theory to 
time-dependent transport in nano-electronic devices 
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In this work a practical scheme is developed for the first-principles study of time-dependent 
quantum transport. The basic idea is to combine the transport master-equation with the 
well-known time-dependent density functional theory. The key ingredients of this paper include: 
(i) the partitioning-free initial condition and the consideration of the time-dependent bias voltages 
which base our treatment on the Runge-Gross existence theorem; (ii) the non-Markovian master 
equation for the reduced (many-body) central system (i.e. the device); and (iii) the construction 
of Kohn-Sham master equation for the reduced single-particle density matrix, where a number of 
auxiliary functions are introduced and their equations of motion (EOM) are established based on 
the technique of spectral decomposition. As a result, starting with a well-defined initial state, the 
time-dependent transport current can be calculated simultaneously along the propagation of the 
Kohn-Sham master equation and the EOM of the auxiliary functions. 
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I. INTRODUCTION 

Quantum transport through nanostructures (e.g. the 
semiconductor quantum dots or organic molecules) would 
play essential role in nano-electronic devices. A rigorous 
treatment should handle the effect of the electronic struc- 
ture of the central device as well as the effect of the inter- 
face to the external contact. This calls for combining the 
theory of quantum transport with the first-principles cal- 
culation of electronic structures. In recent years, consid- 
erable efforts have been focused on the density-functional 
theory (DFT) based simulations on such transport de- 
vices Q, 1^ 1^' In particular, two types of formalisms 
were involved: one is the Lippmann-Schwinger formal- 
ism by Lang and coworkers j^J; the other is the first- 
principles non-equilibrium Green's function (nGF) tech- 
nique [E IE H 0- However, all of these studies 
were restricted to the steady-state transport under time- 
independent bias voltages. 

On the other hand, time-dependent transport phe- 
nomenon are of interest in various contexts, such as the 
single electron pumps and turnstiles, switching-on tran- 
sient behaviors, and ac response in the applications of 
high-frequency amplifiers or detectors. Moreover, it will 
be desirable if one is able to model the nano-electronic 
circuit elements with resistors, capacitors, and inductors. 
To ascribe appropriate quantum capacitances and induc- 
tances to the nano-circuits, the study of time-dependent 
quantum transport is needed. With this motivation, pio- 
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neering work has been carried out [TolllllIT^ ll^ IT^ ITsIl . 

where the time-dependence enters via self-consistent pa- 
rameters and the entire study was largely based on over- 
simplified model description. 

Obviously it is desirable to extend this kind of study 
to the level of first-principles simulation. To this end, 
a necessary element to be combined is seemingly the 
time-dependent density functional theory (TDDFT) |0 , 
which is a generalization of the well-known (ground state) 
static DFT [13 El- Similar to the static DFT, the fun- 
damental variable in TDDFT is no longer the many-body 
wave function but the density. This time-dependent den- 
sity is determined by solving an auxiliary set of nonin- 
teracting time-dependent Schrodinger equations, say, the 
Kohn-Sham equations. The nontrivial part of the many- 
body interaction is contained in the so-called exchange- 
correlation (xc) potential, for which reasonably good ap- 
proximations exist. Since the foundation of TDDFT, an 
enormous amount of progress has been made and the the- 
ory was applied to a large number of problems in physics 
and chemistry. In particular, even the simplest approx- 
imation to the xc potential, i.e., the so-called adiabatic 
LDA, can yield remarkably good results. Nevertheless, 
despite the wide range of applications, the TDDFT has 
been mostly limited to isolated systems. To our knowl- 
edge, its application to quantum open systems (e.g. in 
quantum transport) just begins very recently. These in- 
teresting efforts include either employing formally the 
nGF formalism 20, 22], or propagating wave-function 
at zero temperature j2l| , or introducing electron-phonon 
scattering to balance the external field in order to reach 
a stationary current |2E 01 • 

In this work, by combining the TDDFT with our re- 
cently constructed quantum master equation formalism 
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to quantum transport , we attempt to de- 

velop an alternative scheme. The master equation ap- 
proach to mesoscopic transport can be dated back to 
the classical rate equation '2^ and its quantum gener- 
alization j3^]. More recent work based on the master 
equatio n app roach can be found, for instance, in Refs. 
I31L Is^ l33l 13^ This approach is of interest not only 
by its conceptual difference from the conventional meso- 
scopic transport theory, e.g., the Landauer-Biittiker the- 
ory or the non-equilibrium Green's function (nGF) tech- 
nique [3^ , but also due to its convenience in appli- 
cations. More importantly, since the master equation 
is time-dependent, it seems thus a natural framework for 
the study of time-dependent transport. This is in fact the 
main motivation for us to combine it with the TDDFT 
to construct a time-dependent transport scheme at the 
level of first-principles. 

In the standard TDDFT, the propagation of the sys- 
tem state is described by the time-dependent Kohn-Sham 
equation, whose solution is used to calculate the (time- 
dependent) density and the effective Hamiltonian of the 
non-interacting Kohn-Sham system. In the context of 
quantum transport, using the same idea of TDDFT, we 
first map the entire system (i.e. the central device plus 
the electrodes) to the non-interacting Kohn-Sham sys- 
tem, then trace out the degrees of freedom of the elec- 
trodes. As a consequence, directly related to the time- 
dependent electron density, we obtain a master equation 
for the reduced density matrix of the device, which is a 
counterpart of the well known Kohn-Sham Schrodinger 
equation. We notice that this treatment can be exact 
in principle. That is, the entire system of the device 
plus electrodes evolves under the influence of bias voltage, 
starting from a well-defined initial state. Then, the re- 
quirement of the Runge-Gross existence theorem is satis- 
fied. Of course, as any other (time-dependent) transport 
theories, proper approximate consideration is needed for 
the biased electrodes, which enables us to eliminate the 
degrees of freedom of electrodes and obtain a master 
equation for the reduced state of the central device. 

The remainder of the paper is organized as follows. In 
next section we first specify the transport setup and some 
physical considerations, then outline the main result of 
the master equation approach to transport. In this work, 
instead of the Markovian prescription |25L l2a . 1271 123 | , we 
would like to adopt the scheme of non-Markovian mas- 
ter equation for reasons to be specified later in the main 
text. For completeness, a brief derivation for the non- 
Markovian "n" -resolved master equation is provided in 
Appendix A. In Sec. Ill we combine the TDDFT scheme 
with the multi-particle-state master equation obtained in 
Sec. II. By introducing the single-particle density matrix 
and a number of auxiliary functions associated with the 
spectral decomposition, we establish the central result of 
this work, say, the Kohn-Sham master equation for the 
reduced single-particle density matrix, and the propagat- 
ing equations for the auxiliary functions. These quanti- 
ties would suffice the calculation of the transport current. 



Also, a brief description for the technique of spectral de- 
composition will be presented in Appendix B. In Sec. IV 
we describe how the effect of inelastic electron-phonon 
scattering in the device can be conveniently included into 
the Kohn-Sham master equation. Finally, summary and 
discussions are given in Sec. V. 

II. TRANSPORT MASTER EQUATION 

In this section we shall first present the transport 
model together with some physical considerations on it, 
then outline the main result of the transport master equa- 
tion in many-particle-state Hilbert space, and put its 
derivation in Appendix A. All of these will form the ba- 
sis of the TDDFT scheme to quantum transport of next 
section. 



A. Model Consideration 

Conventionally, mesoscopic transport setup can be de- 
scribed by the transfer Hamiltonian 

H = Hcial,a^)+ ^ ^ea.^.kdl^kdaf.k 

a.—L,R fik 

+ H.C.). (1) 

a—L,R 

Here He is the device Hamiltonian of the central region, 
which can be rather general, e.g., including electron- 
electron and electron-phonon interactions, and may be 
extended to contain a few interface atoms of electrodes 
in the context of such as molecular transport devices. 
The second and third terms describe, respectively, the 
(left and right) electrodes and the tunneling between 
them and the central device. Note that in the above 
Hamiltonian, as widely adopted in transport studies, the 
electrode electrons are treated as noninteracting after 
absorbing the interactions into self-consistent potential. 
This noninteracting feature of the electrodes is in fact 
the basis of the Landauer-type theory, which is the cor- 
nerstone of quantum transport. Also, the electrode elec- 
trons are described in Bloch-state representation. This 
choice is largely for the convenience of formulation. To 
make contact with the first-principles calculation, the lo- 
cal Wannier-state representation will be more appropri- 
ate. As will be shown in Sec. Ill, the conversion between 
them can be easily implemented via the link of self-energy 
matrix. 

Since our major concern is the time-dependent and 
transient behaviors, care should be made for the choice 
of initial condition. In practice, two types of initial states 
were adopted: (i) In the nGF approach to quantum trans- 
port, Caroli et al. considered the two leads as iso- 
lated subsystems with different chemical potentials in the 
remote past. The current starts to flow through the sys- 
tem once the contacts between the device and the leads 
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have been established. This treatment of partitioning is 
seemingly a bit fictitious, since in a real experiment the 
whole system is in thermodynamic equilibrium before an 
external bias is applied deep inside the electrodes. Later 
on this scheme was adopted by Meir and Wingreen et 
al. to obtain the steady-state current through an inter- 
acting central device jSg, and also to time-dependent 
phenomena [s^ . (ii) Conceptually differing from the one 
by Caroli et al., an alternative scheme was suggested by 
Cini in which the central device is contacted and 
in thermodynamic equilibrium before an external time- 
dependent disturbance (i.e. the bias voltage) is switched 
on. This type of initial condition was recently applied by 
Stefanucci ei a? E3,l2Sl23. 



In this work, we involve the Cini's initial condition as 
follows. Initially, before the bias voltage is switched on, 
we let the central device be contacted with the leads and 
in thermodynamic equilibrium. As the bias is switched 
on, the external potential and the disturbance caused 
by the device region are screened deep inside the elec- 
trodes, thus the density inside the electrodes approaches 
the equilibrium bulk-one. This leads to an enormous sim- 
plification since the initial-state self-consistency is not 
disturbed far away from the device. This picture holds 
to be valid provided that the driving frequency is smaller 
than the plasma frequency, which is tens of THz in typical 
doped semiconductors. The bias is established entirely 
across the device by charge accumulation and depletion 
near the electrode-device interfaces. The formation of 
these charge layers then causes a rigid shift of the conduc- 
tion band of the electrodes, but remains the population 
unchanged as the initial one. As an equivalent and more 
convenient description, the above consideration is typi- 
cally replaced by the statement that the reservoir elec- 
trons are always in local thermal equilibrium character- 
ized by the chemical potentials fJ-L/nit) which define the 
bias voltage V{t) via eV{t) — fi^it) — fJ-Rit). In this way, 
the initial equilibrium state is violated, and the system 
starts evolving under the driving of the bias voltage, lead- 
ing to a time-dependent transport current. In practice, 
there exist two types of bias voltages: One is the slow 
time-dependent modulation of the bias voltage, which is 
appropriate for the study of ac response; another is the 
one suddenly switched-on, i.e., eVit) ~ (^l — A*i?)0(^): 
which corresponds to the initial setup for studying the 
transient behaviors. 

Notice that in the above treatment, we explicitly keep 
the chemical potentials (Fermi levels) of the electrodes 
at instant values determined by the time-dependent bias 
voltage, and let the electrode reservoirs always in local 
thermal equilibrium. Physically, this treatment properly 
take into account the closed nature of the transport cir- 
cuit and the rapid relaxation in the electrode reservoirs. 
In Cini's treatment, two spatially homogeneous and time- 
dependent electrical potentials are added to the electrode 
Hamiltonians. The use of the initial occupation in the 
electrodes then corresponds to the initial rigid shift of 
the Fermi levels. However, in the subsequent evolution it 



seems unclear how the rapid relaxation in the electrode 
reservoirs is included, and how the (approximate) local 
thermal equilibrium of the reservoirs is guaranteed. 



B. "n" -Resolved Master Equation and Transport 
Current 

The key idea of the quantum master equation ap- 
proach to transport is to view the device as an open 
dissipative system, and the electrodes as its environ- 
ment. Following the standard treatment of quantum 
dissipation theory, we introduce the reservoir operators 
= Safe taf^kdaf^k = fLf^ + /i?,^- Accordingly, the tun- 
neling Hamiltonian H' reads H' = {o-jiPfj, + H.c.) . 
Let us then consider the entire system evolution start- 
ing from such initial state as discussed above, i.e., the 
central device in a state formed under zero bias in the 
presence of device- electrode coupling, and the electrodes 
in local thermal equilibrium states defined by the initial 
chemical potentials. In the later on evolution, we treat 
H' as perturbation. Note that this does not mean the 
partitioning of the device from the electrodes before the 
bias voltage is switched on. 

Treating H' perturbatively under the second-order 
Born approximation ,28. ^ ..4L .42] . and performing the "n"- 
resolved electrode states average as shown in Appendix 
A, we obtain 



/l(-) „t _ „t Ai+) 



'<i.)HX-«MitVi„^,]+H.c.}, (2) 



where Aj^fj ^ Ea=L,i? ^ith 



^l7.,„pW - E / dt'Ci-l{t,t')g{t,t')[a,p^-Ht')] 
^ Jo 

^ptLpW = E f dt'ciXl{t\t)g{t,t')[p^-\t')aA. 



(3) 



The bath correlation functions read CaiX{t,t') = 

{U{t)fL{t')). and ciXi{t',t) = {fUt')U^.it)). Note 
that here, differing from Ref. we have adopted the 
non-Markovian (i.e. time non-local) form of master equa- 
tion. 

With the knowledge of ^("^(t), one is readily able to 
compute the various transport properties, such as the 
transport current and noise spectrum [2^ |2^ |23, H^- 
In this work, we focus on the calculation of transport 
current. Based on the probability distribution func- 
tion P{n,t) = Tr[p(")(t)], inserting Eq. Q into I{t) = 
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e J2n iT-Pi'iT'i t) gives rise to 



I{t) = 2e ^ ReTr 



(4) 



Here, ^pfl[((0 are the same as defined by Eq. Q with the 
replacement of p^") by the unconditional density matrix 
p = p^"-*. Summing up Eq. (0) over "n", we obtain 

p = -^Cp - ^ { [at , (t) (t)] + H.c. } , (5) 



where = J2a=L,R^pai(t)- In principle, to carry 

out ApQ)i(<) in this non-Markovian scheme, one should 
know p{t') for all the time t' G [0,t), which is in fact 
the non-Markovian nature (i.e. the time non-locality or 
memory effect). However, as will be shown in the follow- 
ing, based on the technique of spectral decomposition, 
we are going to estabhsh the EOM of Alalit). Then, 
the combined propagation of ApQ)i(i) and p{t) will give 
the full time-dependent solution for the transport prob- 
lem under study. Note that these coupled equations will 
be time local, their propagation is thus quite straightfor- 
ward. In contrast, as shown in Ref. 26, if we start with 
a Markovian scheme for the many-particle-state master 
equation, the time-local scheme for the reduced single- 
particle density matrices is seemingly impossible to be 
constructed, and the time nonlocal feature will make the 
practical propagation quite difficult. 



III. TDDFT SCHEME 

The transport master equation constructed in the pre- 
vious section is defined in many-particle Hilbert space. 
This will restrict its applicability to few-states systems, 
because the huge dimensions of many-particle Hilbert 
space for large-scale systems would make the problem 
intractable. In this section, in the spirit of TDDFT we 
recast the many-particle interacting system to a Kohn- 
Sham noninteracting system, and establish the corre- 
sponding transport master equation, which is defined in 
the single-particle-state Hilbert space with dimensions 
greatly reduced. 



A. General Consideration 

Taking the entire system of electrodes plus device into 
account and starting with the well-defined initial state 
as described in Sec. 11(A), the Runge-Gross theorem im- 
plies a one-to-one correspondence between the electron 
density and the potential function. Then, the electron 
density distribution uniquely determines all the (trans- 
port) properties of the system. The significant role of 
the Runge-Gross existence theorem is to guarantee the 



construction of a proper noninteracting Kohn-Sham sys- 
tem. 

For the electrodes, the Kohn-Sham mapping to non- 
interacting system is relatively simple. In practice, the 
interactions in the electrodes are absorbed into a self- 
consistent potential, and the electrons in the electrodes 
are treated as noninteracting. As a matter of fact, the 
noninteracting feature of the electrodes is the basis of the 
Landauer-type transport theory, which is the cornerstone 
of quantum transport. 

The part of the system that we are going to treat care- 
fully is the (extended) device. Its Kohn-Sham counter- 
part can be constructed as follows. In principle, if we 
knew the many-particle density operator p{t) of the de- 
vice, we could introduce the reduced single-particle (RSP) 
density matrix, cr^^{t) = Tr[ata^p(t)]. However, rather 
than solving /o(t), we should calculate afj_^{t) directly 
from its equation of motion, which is the central gaol of 
this work. Within the TDDFT framework 16], we map 
the device interacting Hamiltonian to the non-interacting 
Kohn-Sham one as 



hrnnit) = hl^^{t) + v'^Jt) + ^ a,j{t)Vm 



(6) 



In first-principles calculation the state basis 
is usually chosen as the local atomic orbitals, 
{(/)m(r),m = 1,2,---}. Here hP{t) is the non- 
interacting Hamiltonian which can be in general time- 
dependent; Vmnij is the two-electron Coulomb integral, 
Vmmj = /dr/drV;,(r)0„(r)^0*(r')0j(r'); and 
v^nit) = /dr</.;,(r)^-[n](r,i)<^„(r), with v-'=[n]{r,t) 
the exchange-correlation potential, which is defined 
by the functional derivative of the the exchange- 
correlation functional A^'^. In practice, especially 
in the time-dependent case, the unknown functional 
A^'^ can be approximated by the energy functional 
E'^'^, obtained in the Kohn-Sham theory and fur- 
ther with the local density approximation (LDA). 
Note that the density function n(r, t) appearing in 
the Kohn-Sham Hamiltonian is related to the RSP 
density matrix via n{r,t) = '^"l('')^™"(*)C(I■)■ 
Thus, the entire Kohn-Sham system is described by 
the noninteracting electrodes, the device Hamiltonian 
H{t) = Y,mnhmnit)alnan, and the coupling between 
them. Starting with this entire Kohn-Sham Hamiltonian 
and repeating the derivation in Sec. II (B) will lead 
to the same formal result presented there. The only 
difference is that now the Kohn-Sham Hamiltonian H{t) 
is noninteracting, which enables us to establish the EOM 
for (T^^(t). 



B. Kohn-Sham Master Equation 

In this subsection, for the Kohn-Sham system, we re- 
cast the master equation Eq. lO into the EOM of cr^^it) 
and a set of auxiliary functions, and re-express the trans- 
port current in terms of them. For brevity, we tenta- 
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tively restrict our derivation to the special case of con- 
stant bias voltage. Generalization to time-dependent 
voltages is straightforward, and will be discussed in Sec. 
Ill (D). For constant bias voltage, the correlation func- 
tions c'"ai}u{t,t') and C^lt{t',t) only depend on the dif- 
ference of times. By employing the technique of spectral 
decomposition as described in Appendix B, they can be 
formally decomposed in the sum of exponential functions 

CUM = Y^xU^e^U'i^-^'\ 

k 
k 

where the parameters A^^-**^ and ^^^^'^ are uniquely de- 
termined by the spectral decomposition. Furthermore, in 
addition to a^f^{t), we introduce 

a,^{t,t') - Tr{atg(t,t')[a,p(t')]}, (8a) 
a,^{t,t') = TY{alg{t,t'Mt')a,]}, (8b) 

and the auxiliary functions 

Ai-JX'it) = /'dt'AL;l'^e-^.'^^(-*')a.,,(M'), (9a) 
Jo 

AiV'At) - /'d^'AW>--(*-*')^..'(i,0. (9b) 
Jo 

Simple algebra leads to the EOM of (j{t, t') and a{t, t'): 

dtaut,{t,t') = t[a{t,t')h{t)]^^, (10a) 
dta,t.{t,t') = i[a{t,t')h{t)]^^. (10b) 

Accordingly, the EOM of A'^^J^^, (t) can be straightfor- 
wardly obtained: 

dtAl~JX,{t) = Xi-i'a,,'{t) 

+7i72'4;l';'(i) +*E^U™W'^™m'(0, (11a) 

dtAi%{t) = xMk[s,^,^a,^,(t)] 
+7itl.'=4tl'.'W +*E^"U™W'^™m'(0- (lib) 

m 

Note that the RSP density matrix a^)j,{t) appears in the 
EOM of these auxiliary functions. To close the EOM 
for the transport problem, we need to derive the EOM 
of (J^f_,{t). Using the identity Tr{at a^, X)^' [a|i', A^-]} = 
Tr[at Ajy], we easily obtain 

<7{t) = -i[h{t),a{t)] - E t*^" W + H-c] , (12) 

a 

where the matrix AIa{t) is defined via its elements as 



Propagating Eqs. Hll|l and (|12|l is numerically straightfor- 
ward, and the most difficult task arising from the mem- 
ory kernel in usual non-Markovian dynamics is avoided. 
Finally, the transport current is simply related to the 
auxiliary functions as 



I{t) 



2e E R-e [. 



A 



(-)fe 

Rfivfj, 



(i)-<l',W 



(13) 



which is output automatically along the time propagation 
of the above EOM. 



C. Spectral Density Function 

The treatment in the previous subsection is 
based on the spectral decomposition of r^^(e) = 
^QA'fe^aiyfe^(e — efc), as shown in Appendix B. Since 
the spectral decomposition is rooted in a technique of 
numerical fit, it is thus in principle suited for arbitrary 
shape of the spectral density function r^j^(e). In par- 
ticular, it is not limited by the conventional wide-band 
approximation. This advantage makes us be able to 
incorporate the electronic structure of the electrodes 
which is obtained from other sophisticated methods. 
One possible way of calculating r^^,(e) is based on a 
semi-empirical tight-binding model for the electrodes, 
and using the surface Green's function technique. Some 
details of this method can be found in Ref. Here we 
would like only to outline the key idea for completeness. 

Formally, in matrix form, the spectral density function 
is related with the self- energy matrix via 



(14) 



Further, the self-energy matrix can be calculated by 

T.^tgt\ (15) 

where "t" is the coupling matrix between the "edge 
atoms" of the electrode and the (extended) device, and 
g is the surface Green's function of the electrode. Recur- 
sively, g can be obtained via the following Dyson equation 
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9o 



9o 



igi\ 



(16) 



where "i" is the coupling matrix between nearest- 
neighbor layers in the electrode. 

In practice, care should be paid to the dimensions of 
the matrices, i.e., there are two types of orbital labels: 
one is restricted to the device edge-orbitals; and another 
is over all the device orbitals. From the definition of each 
matrix, the size of its every dimension can be identified 
accordingly. 



D. Time-Dependent Voltage 

The derivation in Sec. Ill (B) corresponds to time- 
independent bias voltage. That is, after a constant volt- 
age is switched on, the chemical potentials in electrodes 
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remain unchanged. This imphes that the correlation 
functions ci^^(t,t') are of time-translational invariance 
during the later evolution. The description under this 
assumption can be applied to the study of transient be- 
haviors [3^ . 

Another important time-dependent setup is applying 
modification on the bias voltage, such as in the study of 
ac response. As have discussed previously, under proper 
conditions, the effect of time-dependent voltage can be 
approximately described by rigid shifts of the conduction 
bands of the electrodes, i.e., eka{t) = ^ka + ^a{t), and 
keeping the occupation on each state unchanged. Ac- 
cordingly, the correlation functions read 



Formally, the many- electron- state master equation for 
the device can be expressed as 



Cai^fi it , ^) 



Ca^iy (t i ) 
Coif^ {t i ) 



-I //, titiA„(ti) 



(17) 



where C^u{t ~ t') are the counterparts in the absence of 
time-dependent voltage shift (i.e. ^aii) — 0), and can be 
obtained by using the surface Green's function technique 
as described in the previous subsection. 

Replacing the correlation functions in Sec. Ill (B) with 
Eq. H17|l . all the equations obtained there can be for- 
mally re-derived, except with only a minor modifica- 
tion on the second terms of the l.h.s. of Eq. i.e.. 



'•^■''^ — i^aif)- The advantage of the proposed 



TDDFT scheme in Sec. Ill (B) is thus prominent, since in 
other conventional treatments the time-dependent trans- 
port is usually regarded more difficult than its stationary 
counterpart, owing to the lack of time-translational in- 
variance. However, in our scheme, no extra efforts are 
needed for time-dependent voltage in propagating the 
EOM of the RSP density matix a^,j{t) and the auxiliary 

functions vl^^lf (i). 



IV. ELECTRON-PHONON INTERACTION 

In this section we consider the issue of inelastic scat- 
tering in the device. Still within the TDDFT framework, 
we assume the Kohn-Sham subsystem of the device being 
coupled to a phonon bath (-ffph = 'Y^q^qb\bq)- In gen- 
eral, the electron phonon interaction Hamiltonian has the 
form 

q mn 

^Y.^wji + wy,^), (18) 



qK 



where for brevity we introduce = Wmn = aJnOn, and 
fqK is defined accordingly. Here and in the remainder of 
this section, we adopt the electronic eigenstate basis of 
the device Kohn-Sham Hamiltonian before the bias volt- 
age is applied. This choice has the advantage of making 
the electronic state transition due to phonon scattering 
clearly defined. 



p = -iCp - UcP - TZphP- 



(19) 



In this equation, the term TZcP describes the electrode 
effect on the device, and TZphp stems from the effect of 
electron-phonon interaction. In the previous sections, we 
have focused on the term TZeP, by performing a non- 
Markovian treatment at the level of second-order Born 
approximation. For TZpiip, one can in principle perform 
the same treatment, by also applying the spectral decom- 
position technique for the phonon bath. Nevertheless, for 
simplicity we would like to treat the electron-phonon in- 
teraction under the Markovian approximation as usual. 
Following Ref. we have 

T^phP = I E {[^«' ^^^-^ - -^^^^'l + H-C-} • (20) 



At the transition energy lj^ = |em ~ Enl, the operators 
W!^^^ read w!^^^ = vi^^W^, where fi^^ are defined by 
fi+'> = 5(c^«)|7«|2(n^^ + 1) and vi'^ = g{u,)\j,\^fiu,^, 
respectively. Here g{uj) is the density of states of the 
phonon modes, and is the corresponding phonon occu- 
pation number. Also, we would like to mention that to ar- 
rive at Eq. (|20|l . we have inserted the device Kohn-Sham 
Hamiltonian at the initial equilibrium state, but not the 
time-dependent one associated with the later evolution, 
into the dissipation terms. This is the well-know approx- 
imation in studying dissipative systems under external- 
field driving. This treatment reduces Eq. (|20|l to the 
Lindblad form. 

To incorporate the effect of electron-phonon interac- 
tion into the Kohn-Sham master equation (|12|l . we need 
to recast TZphP to a RSP density matrix form. Simple 
algebra gives 



n n 

(21) 

In the derivation of this result, the Wick- type factor- 
ization such as (aj, a„Wrii/) ~ av^ann is assumed. We 
would like to note that Eq. (|21|l coincides precisely with 
the central result derived by Gebauer and Car in an alter- 
native transport approach |24j . After simple basis trans- 
formation (i.e. from eigenstates to local atomic orbitals 
) , inserting Eq. (|21|l into the Kohn-Sham master equation 
(|12|l leads to an elegant scheme which can also account 
for the electron-phonon scattering in the time-dependent 
transport process. We stress that this is another signifi- 
cant advantage of the proposed transport approach. 
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V. CONCLUSION AND DISCUSSION 

To summarize, based on the quantum master equation 
approach, we have constructed a practical scheme for the 
first-principles study of time-dependent quantum trans- 
port. The basic idea is to combine the transport master- 
equation with the well-known time-dependent density 
functional theory. By utilizing the partitioning-free ini- 
tial condition, the scheme is reliably based on the Runge- 
Gross theorem. Then, with the help of spectral decom- 
position, a closed set of equations for the RSP density 
matrices are obtained. Time propagation of this set of 
equations will directly output the time-dependent trans- 
port current. 

It is of interest to note that it is the non-Markovian 
(but not the Markovian) master equation that makes 
us be able to construct the closed set of equations for 
the RSP density matrices, based on the spectral decom- 
position technique. In Ref. 129 . without using the spec- 
tral decomposition and in the framework of Markovian 
master equation, we established an alternative form of 
TDDFT based master equation for the RSP density ma- 
trix. However, in that scheme, the numerical propaga- 
tion seems inefficient, because of the non-local effect of 
time integration. On the contrary, the numerical imple- 
mentation of the present scheme should be efficient and 
straightforward. Moreover, it will be very flexible for the 
parametrization of highly structured spectral densities, 
which makes the description far beyond the Markovian 
approximation. Systematic applications and numerical 
implementations of the proposed scheme are in progress 
and will appear in the forthcoming publications. 

Finally, we remark that the major approximation in- 
volved in the master equation is the second-order Born 
approximation for the tunnel Hamiltonian. This is a 
standard and well-justified approximation, which makes 
the resultant master equation good enough in a large 
number of dissipative systems (e.g. in quantum optics). 
Also, our recent work clearly showed its satisfactory ac- 
curacy in quantum transport [2^ [SJ, H^. In this 
context, we may roughly claim that its accuracy is at 
least at the level of sequential transport. Noticeably, 
in the first-principles study, other numerical errors will 
completely cover up the inaccuracy of sequential trans- 
port. Therefore, the proposed TDDFT master equation 
scheme should be an attracting theoretical tool for the 
first-principles study of time-dependent transport. 

APPENDIX A: "n"-RESOLVED MASTER 
EQUATION 

In this appendix we present the derivation of the 
non-Markovian form of the "n" -resolved master equa- 
tion. In almost the same spirit of Ref. let us regard 
H' — {o'ltPfj. + H.c.) as a coupling of the system of 
interest to a dissipative environment. Treating H' as per- 
turbation and up to the second order, a formal equation 



for the reduced density matrix is obtained as |4l| 

p(t) = -Z£p(t) - f dT{C'{t)g{t,T)C'{T))p{T). (Al) 

Jo 

Here the Liouvillian superoperators are defined as 

£(...) ^ [i?s,(---)], ^'(•■•) = [H'. {■■■)], and 
g{t,T){---) = G{t,T){---)G'<{t,T) with G{t,T) the usual 
propagator (Green's function) associated with the cen- 
tral system (device) Hamiltonian Hs- The reduced den- 
sity matrix p{t) = TrB[/9T(t)], and ((• • • )) = TrB[(- • • )pb] 
with pb the density matrix of the electrode reservoirs. 

The trace in Eq. (|A1|) is over all the electrode degrees 
of freedom, leading thus to the equation of motion of 
the unconditional reduced density matrix of the central 
system. However, more information is to be contained if 
we keep track of the record of electron numbers arrived 
at the collector (right electrode). We therefore classify 
the Hilbert space of the electrodes as follows. First, we 
define the subspace in the absence of electron arrived at 
the collector as "5'^°)", which is spanned by the product 
of all many-particle states of the two isolated reservoirs, 
formally denoted as _B(°) = span{|^'i) (gi {"^r)}. Then, 
we introduce the Hilbert subspace "S'"^" ( n = 1, 2, • • ■ ), 
corresponding to "n" electrons arrived at the collector. 
The entire Hilbert space of the two electrodes is B = 

With the above classification of the reservoir states, 
the average over states in the entire Hilbert space "B" 
in Eq. Q is replaced with states in the subspace 
leading to a conditional master equation 

p(")(t) = -iCp^^Ht)- f dTTYst„,[c'{t)g{t,T) 

Jo 

xC'{t)pt{t)]. (A2) 

Here p^"^\t) = Tr^c,!) [pt(0]; which is the reduced den- 
sity matrix of the system conditioned by the number of 
electrons arrived at the collector until time t. Now we 
transform the Liouvillian operator product in Eq. (|A2I) 
into the conventional Hilbert form: 

C'{t)g{t,T)C'[T)pT{T) 
= [H'{t)G(t,T)H'{T)pT{T)GHt,T) 

~G{t, t)H'{t)pt{t)G^ (t, T)H'{t)] + H.c. 
= [/ - //] + H.c. (A3) 

To proceed, two physical considerations are further taken 
into account as follows: (i) Instead of the conven- 
tional Born approximation for the entire density ma- 
trix pt{t) — p{t) ® Pb, we propose the ansatz ptIt) — 

J2n p'"'' (''') ® p'b^ ' where is the density operator of 
the electron reservoirs associated with n-electrons arrived 
at the collector. With this ansatz for the density opera- 
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tor, tracing over the subspace "_B(")" yields 

x[a^G(t,r)atp(")(r)Gt(t,r)] 

+T:TB[F,{t)Fl{T)p^;^^] 

x[atG(i,r)a,p(")(r)Gt(t,r)]} 

(A4a) 

Tr5(„,[//] = E{TrB[/i.(r)p(j"VL^W] 

+TrB[/L.(r)pi"Vi^W] 
x[G(t,r)atp(")(r)Gt(t,r)a^] 

+TrB[/L(T)p^'"'VHM(i)] 
x[G(t,r)a,p("-i)(r)Gt(t,r)at] 

x[G(t,r)4p("+i)(r)Gt(t,T)a^]}. 

(A4b) 

Here we have utilized the orthogonality between states in 
different subspaces, which in fact leads to the term selec- 
tion from the entire density operator p^. (ii) Due to the 
closed nature of the transport circuit, the extra electrons 
arrived at the collector will flow back into the emitter 
(left reservoir) via the external circuit. Also, the rapid 
relaxation processes in the reservoirs will quickly bring 
the reservoirs to the local thermal equilibrium state de- 
termined by the chemical potentials. As a consequence, 
after the procedure (i.e. the state selection) done in 
Eq. HA4|I . the electron reservoir density matrices and 

p^^^^ should be replaced by i.e., the local thermal 
equilibrium reservoir state. Then, the non-Markovian 
"n" -resolved master equation Eq. Q is obtained. 

APPENDIX B: SPECTRAL DECOMPOSITION 

In this appendix we show how to express the cor- 
relation functions C^t{t,t') — {fafj,it)fl^{t')) and 

Git^(t',t) — {fi^{t')fafi{t)) in terms of the sum of ex- 
ponential functions, via the technique of spectral decom- 
position. For constant bias voltage, let us re-express, for 
instance, C^li{t',t) as 

k 

- / ^r^,(e)n„(6)e-(*-*'), (Bl) 

where r^^(e) = 27r X^a, iLfc^aMfc^(e - Cfc) is the spectral 
density function of the a-th electrode, which can be ob- 



tained by the surface Green's function technique as de- 
scribed in Sec. III(C). 

To express Git^(i',t) as a sum of exponential func- 
tions, let us parameterize the spectral density as follows 

c.w = E(^_^,j;; (rL.)^' ^^^^ 

By employing the theorem of residues, the integral result 
of Eq. (|BT1) reads 

CLtlit\t) = ^|^n„(f2L.)e-^^^^''(*-*') 
k=i "•'t^ 

OO 

+ ^Er".(^n)e"""(*-*'\ (B3) 

where fl^^^ = - iV'^^^, and e„ = ep - i{2n + l)/f3. 

Similarly, the correlation function Gi^i(t,t') can be de- 
composed. Then, these two functions can be formally 
re-expressed as Eq. 0. 

About the spectral decomposition, a few remarks are 
made as follows: (i) This spectral decomposition tech- 
nique has been successfully applied to the non-Markovian 
dissipative systems E4| . and to the quantum transport 
through molecular wires very recently |34| . The most 
prominent advantage of this technique is its flexibility. 
That is, it is extremely well suited for the parametriza- 
tion of highly structured spectral densities, leading to 
long and oscillatory correlation functions, thus making 
the description far beyond the Markovian approxima- 
tion. In the context of quantum transport, this decom- 
position technique avoids the usual wide-band approxi- 
mation, and can handle arbitrary band structures which 
may be obtained even by the first-principles calculation, 
(ii) In reality, we might find different sets of parame- 
ters which can approximate a given spectral density to 
the same degree of accuracy. Nevertheless, this will not 
lead to a different dynamical behavior, since the dynam- 
ics is determined by the spectral density itself, (iii) Not 
all types of spectral densities need the numerical decom- 
position of the spectral densities. For instance, for the 
Drude spectral density J{uj) — riuj/[l + {u!/ud)^], or other 
spectral densities with simple poles, there is no approxi- 
mation other than the finite number of Matsubara terms 
involved. In practice, the number of fitting terms de- 
pends on the shape of spectral density. For example, 
for the Ohmic spectral density with exponential cutoff 
J{uj) = r\LLie~^l^'' ^ only three terms are accurate enough 
to parametrize the spectral density. However, for another 
spectral density of the form J(w) = 7/u;^/(2ci;^)e~'^/'^'=, 
nine terms were necessary for an accurate fit j44j . (iv) 
For finite temperatures, the contributions of the high 
Matsubara frequencies will be small, so that for practical 
purposes the summation over the Matsubara frequencies 
can be truncated at some value. Hence, this decomposi- 
tion will be most profitable for high temperatures, when 
only a few Matsubara frequencies contribute. For very 
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low temperatures, we may alternatively parametrize the 
combined spectral density, say, r"^(e) = T'^^{e)na{e), in 
terms of the Lorentzian spectral function Eq. ljB2p . 
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